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A discontinuous Galerkin (DG) method is formulated, implemented, and tested for sim- 
ulation of compressible turbulent flows. The method is applied to turbulent channel 
flow at low Reynolds number, where it is found to successfully predict low-order statis- 
tics with fewer degrees of freedom than traditional numerical methods. This reduction 
is achieved by utilizing local hp-refinement such that the computational grid is refined 
simultaneously in all three spatial coordinates with decreasing distance from the wall. 
Another advantage of DG is that Dirichlet boundary conditions can be enforced weakly 
through integrals of the numerical fluxes. Both for a model advection-diffusion problem 
and for turbulent channel flow, weak enforcement of wall boundaries is found to improve 
results at low resolution. Such weak boundary conditions may play a pivotal role in wall 
modeling for large-eddy simulation. 


1. Introduction 

In this paper we formulate, implement, and apply a discontinuous Galerkin (DG) 
method for the simulation of compressible turbulent flows. Discontinuous Galerkin can 
be thought of as a hybrid of finite- volume and finite-element methods that has a number 
of potential advantages including: high-order accuracy on unstructured meshes, local hp- 
refinement, weak imposition of boundary conditions, local conservation, and orthogonal 
hierarchical bases that support multiscale turbulence modeling (Hughes et al. 2000; Collis 
2001, 2002). The interested reader should consult the review of Cockburn (1999) and 
Cockburn et al. (2000) for a recent update on the status of discontinuous Galerkin. Since 
the DG method is ideally suited for hyperbolic or nearly hyperbolic systems, we believe 
that DG may be a particularly attractive method for high-Reynolds-number compressible 
turbulent flows in complex geometries. This paper takes a first step in applying DG 
to turbulent flows by considering low-Reynolds-number DNS of compressible turbulent 
channel flow. We note, before proceeding, that there is considerable ongoing research on 
DG methods (see Cockburn et al. 2000) and we have greatly benefited from the work of 
Cockburn and co-workers, Karniadakis and co-workers, and Bassi and Rebay. 


2. Formulation 

Consider the compressible Navier-Stokes equations in strong form 

U, t + Fi,i - Fli = S in D, (2.1a) 

U(x,0) = Uo{x) at t = 0, (2-16) 

where U = {p, pu, pe} T is the vector of conserved variables, p is the fluid density, u is the 
fluid velocity vector, and e is the total energy per unit mass. The inviscid and viscous 
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== fi i -f- fi 2 



Figure 1 . Schematic of DGM discretization 

flux vectors in the zth coordinate direction are Fi(U) and and S is a source 

term, including body forces in the momentum equations and a heat source in the energy 
equation. Equation (2.1 a) is solved subject to appropriate boundary conditions, which 
must be specified for each problem of interest; a state equation, such as the ideal gas 
equation; and constitutive laws that define fluid properties such as viscosity and thermal 
conductivity as functions of the conserved variables. Due to space limitations, we do not 
explicitly define the flux vectors, state equation, or constitutive relations, but instead 
refer the reader to standard texts such as Hirsch (1988). 

The fixed spatial domain for the problem is denoted by fi, which is an open, connected, 
bounded subset of R d , d = 2 or 3, with boundary 3fi. Let Vh be a partition of the domain 
fi into N subdomains fi e where 

l v 

n=\Jn e and =0 for e^f. (2.2) 

e— 1 

Starting from the strong form of the compressible Navier-Stokes equations (2.1a), we 
consider a single subdomain, fi c , multiply by a weighting function W which is continuous 
in fi e , and integrate the flux terms by parts 

f ( W T U,t + W^iF* - Fi )) dx+ f W T ( F n - F v n ) ds= j W T S ds (2.3) 
n. di. n. 

where F n = If the solution were assumed to be continuous and this equation 

were summed over all the elements in 'Ph ? then all the flux terms would telescope to 
the boundary c?fi and we would obtain the standard continuous Galerkin form of the 
compressible Navier-Stokes equations. However, in discontinuous Galerkin, one instead 
allows the solution and weighting functions to be discontinuous across element inter- 
faces (see figure 1) and the solutions on each element are coupled using appropriate 
numerical fluxes for both the inviscid flux F n (U) ->■ F n (U~ ,U + ) and the viscous flux, 

F\{U,Uj) -* Fi(U~ ,U~,U + ,Uj). Introducing numerical fluxes and summing over 
all elements yields 

£ / (w T U, t + W T i (F v i -F i )) dx + 

e=l A 

N r / Qe N r 

I W T [K{U-,U + )-F V n (U-,U-,U + ,U+)) = W T Sds (2.4) 

e=1 «=i n e 

where the U + and U~ states are defined in figure 1. For an element edge on the physical 
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boundary dfi, U + = Ubc ♦ Likewise, for inter-element boundaries, comes from the 
neighboring element. Thus, all interface and boundary conditions are set through the 
numerical fluxes. Rewriting (2.4) in a more compact notation, the discontinuous Galerkin 
method is: 

Given U 0 = C7 0 (x), for t £ (0 ,T), find U(x,t) € V{V h ) x H l {0,T) such that l/(x,0) = 
Uo(x) and 

B dg (W , 17) = (W, S) VW € V(Pfc), (2.5) 

where V('P^) is the broken space defined in Baumann & Oden (1999). If V(Vh) is restricted 
to a space of continuous functions, then one recovers the classical continuous Galerkin 
approximation upon using the consistency properties of the numerical fluxes (Cockburn 
1999). 

While there is a wide range of choices for both the inviscid and viscous numerical 
fluxes (see Cockburn (1999) for a thorough review), we have initially chosen to use a 
Lax-Friedrichs method for the Euler flux 

Fn(U~,U + ) = | {F n (U~) + F n (U + )) + A m {U- - U + ) (2.6) 

where A m is the maximum, in absolute value, of the eigenvalues of the Euler Jacobian 
A n = dFJdU. 

For the numerical viscous flux, we use the method of Bassi & Rebay (1997). First, a 
“jump savvy” gradient of the state, <jj ~ U j is computed by solving 


f; fv T ajdx = -^ [VjUdx + 52 [ V T U nj ds VV e V{V h ) (2.7) 
e=1 e=1 n. e=1 an. 


for each direction, j , where 

U=l(U-+ U + ) . (2.8) 

The Bassi-Rebay viscous flux is then computed using 


r n (U-,<Tj,U + ,<Tp = \(F v n (U-,'T-) + Fl(U + ,cr+)) . (2.9) 

While this method is known to be only “weakly stable,” (Arnold et al. 2002) we have 
not encountered any difficulties for the problems considered here, and this method has 
been used successfully in the past (Bassi & Rebay 1997). In the future, we will consider 
other, provenly-stable, numerical fluxes for the viscous terms, and the reader is referred 
to Arnold et al. (2002) for an extensive discussion of the advantages and disadvantages 
of a wide range of viscous fluxes for use in discontinuous Galerkin discretizations. 

In setting boundary conditions weakly through the numerical fluxes, one must con- 
struct a state, Ub c , that enforces the appropriate boundary conditions, and Atkins 
(1997) provides a discussion of the important issues involved in selected Ubc ■ For the 
Navier-Stokes calculations reported here, we use the following approach. At far-field 
boundaries Ubc is se t to freestream values. At isothermal wall boundaries, we evaluate 
Ubc separately for the convective and viscous fluxes. Let q\ = (u~n y - v~n x )n y and 
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Q 2 = ( v~n x — u~n y )n x then the reconstructed state at a wall for the convective flux is 


U bc = 


P~ 

P~Qi 

P~Q2 

p~e~ +0.5p~(ql +ql) j 


( 2 . 10 ) 


This state enforces the no-penetration condition which is appropriate for both inviscid 
and viscous calculations. For the viscous flux, the no-slip condition is enforced using 


U bc = 


P 

0 

0 

P“2V(7(7-1)M 2 ) , 


( 2 . 11 ) 


where T w is the prescribed wall temperature, 7 is the ratio of specific heats, and M is 
the reference Mach number. 

By way of summary, the discontinuous Galerkin method is a hybrid of finite-element 
and finite- volume methods, where solutions are continuous within an element but dis- 
continuous across element interfaces, and elements are coupled via numerical fluxes on 
element interfaces. Discontinuous Galerkin has several potential advantages including: 
(1) Spectral accuracy on arbitrary meshes, (2) Local /ip-refinement, (3) Boundary con- 
ditions are imposed weakly through numerical flux, (4) Local conservation allows for 
different fidelity models on neighboring elements, (5) Orthonormal hierarchical basis on 
each element readily supports multiscale turbulence models, and (6) DG works best near 
the hyperbolic limit making it potentially valuable for high Reynolds number turbulence. 
A thorough review of the DG method is available (Cockburn 1999) while a more com- 
plete description of DG for turbulence simulation including a discussion of multi-scale 
turbulence modeling is given in (Collis 2002). 


3. Discretization and implementation 

For every element fi e £ 'Ph we define the finite-dimensional space P Pe (Cl) of polynomials 
of degree < p e defined on a master element Cl. Then 


PpAfte) = {d>\4> = e P p .(fi)} 


(3.1) 


where Jn t is the Jacobian of the transformation of element fl c to the master element 
and 


V P (V h ) = 



c V(V h ) 


(3.2) 


where m is the number of conserved variables; m = 5 in three dimensions. 

Thus, the semi-discrete discontinuous Galerkin method is: Given Uq = U o(x), for 


t £ (0,T), find Uh(x,t ) £ Vp(Vh) x H l (Q,T) such that 


B DG (W h , U h ) = (W h , S) , VW h £ V p (P h ) . (3.3) 

We utilize the family of orthogonal, hierarchical bases formed from tensor products of 
Jacobi polynomials, as described in Karniadakis & Sherwin (1999), which are supported 
in a wide range of elements types in two and three dimensions. For time advancement, 
we currently use the third-order TVD-RK method (Shu 1988; Shu & Osher 1988) 
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Figure 2. Typical parallel speedup for DG implementation on a Pentium IV Beowulf cluster. 

The DG formulation presented above has been implemented using object-oriented pro- 
gramming (OOP) in fully modern ANSI/ISO C++ using the Standard Template Library 
and generic programming concepts. For efficiency, all kernel computations are performed 
using the ATLAS library, and the code runs on a number of operating systems includ- 
ing Linux, Windows, and SGI Irix. The code is implemented as an element library that 
supports all the operations required for discontinuous Galerkin, and we have used this 
library to implement specific solvers for advection-diffusion, Burgers, wave, linearized- 
Euler, Euler, Navier-Stokes equations. Due to the inherent locality in the discontinuous 
Galerkin discretization, parallel implementation is particularly easy and efficient. We use 
the MPI-2 library (including parallel MPI-IO) and parallel efficiency results are shown 
in figure 2 for our Pentium IV Beowulf cluster demonstrating excellent scaling. 

4. Weak boundary conditions 

One of the issues that arises in using discontinuous Galerkin methods is that Dirichlet 
boundary conditions are most naturally enforced weakly through the numerical fluxes. 
While similar “weak” boundary conditions have been used for far-field nonreflecting 
boundary conditions in finite-difference discretizations (see e.g. Poinsot & Lele (1992); 
Thompson (1987)) the use of weak boundary conditions for wall-type boundary condi- 
tions is less common, especially in the flow physics community. In the computational 
mechanics and applied mathematics communities there has been prior work on weak 
enforcement of Dirichlet boundary conditions in the continuous Galerkin method by 
Babuska (1973) and Nitsche (1971) and these methods are related to discontinuous 
Galerkin (Arnold et al. 2002). Likewise, the recent work of Layton (1999) provides the- 
oretical considerations for weakly enforced Dirichlet boundary conditions for the Stokes 
problem that are motivated by observations of improved solution quality compared to 
hard Dirichlet boundary conditions. 

While one can always set “hard” Dirichlet boundary conditions in any discretization 
(including DG), it is interesting to compare the performance of hard boundary conditions 
with weak enforcement through the numerical fluxes as described above. As an example, 
consider the simple steady forced advection-diffusion problem 

U, x = 1 + l 'U,xx (4.1) 

with boundary conditions u(0) = u(l) = 0 and diffusivity v = 0.01. Figure 3 shows re- 
sults computed using a discontinuous Galerkin discretization with two p = 6 elements and 
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X X 

Figure 3. Weak (a) and hard (b) Dirichlet boundary conditions for an advection-diffusion problem 


BC Loo L 2 Hi 
Weak 0.374 0.0198 2.00 
Hard 0.251 0.0850 3.35 

Table 1. Errors in advection diffusion solutions 


both hard and weak enforcement of the Dirichlet boundary conditions. This discretiza- 
tion was intentionally selected to be coarse in order to highlight the differences between 
the two solutions. One clearly sees that oscillations are more pronounced when a hard 
boundary condition is used. Conversely, while oscillations are less in the weak case, the 
boundary condition on the right side (inside the boundary layer) is only approximately 
satisfied; u( 1) = 0.374 instead of zero. Table 1 compares the error in the solution in the 
Loci Li, and Hi norms. Consistent with the graphical results, the solution with weak 
Dirichlet boundary conditions has four times less error in and is also better in H \ . 
However, the solution with weak boundary conditions is slightly worse in L ^ and this is 
directly due to the error in the boundary value. Discarding a small region near x = 1, the 
weak solution is also better in L^. While these results are certainly not conclusive, they 
are indicative of the potential benefit to be gained from weak enforcement of Dirichlet 
boundary conditions that are naturally obtained from a DG discretization. Philosoph- 
ically speaking, one should not enforce boundary conditions any more accurately then 
the error in the interior solution. Doing so tends to over-constrain the interior solution, 
typically leading to oscillations as seen in figure 3(b). By weakly enforcing boundary 
conditions one obtains solutions that still feel the influence of the boundary through the 
numerical fluxes, but in a manner that is consistent with the accuracy of the interior 
solution, leading to improved solutions away from the wall. Given the importance of wall 
boundary conditions for near-wall turbulence, we will pay particular attention to the 
success of the weak boundary condition throughout the following discussion. 


5. Flow over a circular cylinder 

Before applying our DG formulation to a turbulent flow, we begin by considering a 
benchmark problem of both steady and unsteady flow over a circular cylinder. 
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DG 

DG 

FD 

Experiment 

Re 

P 

= 4 

V 

= 6 

6 th ' 

order 




s/d 


s/d 

Cd 

s/d 

c d 

s/d 

c d 

20 

0.96 

2.125 

0.96 

2.124 

0.93 

1.98 

0.9 

2.01 

40 

2.39 

1.589 

2.39 

1.589 

2.36 

1.50 

2.1 

1.48 


Table 2. Drag and separation length for laminar flow over a circular cylinder, with comparison to prior 
computations and experiments. The computational and experimental data are taken from Visbal (1986). 


5.1. Steady flow 

Consider the steady, laminar flow of air past an isothermal circular cylinder kept at 
the freestr eam temperature. The freestream Mach number is M = 0.2 and results are 
reported for two Reynolds numbers: 20 and 40. By considering a series of different domain 
sizes, we eventually selected a domain of fi = [-15, 30] x [-30, 30] as sufficiently large to 
prevent adverse influence on the net drag and length of the separation bubble. A block 
structured mesh using 812 quadrilaterals was generated using a special purpose grid 
generator (Tezduyar 1991) and each quadrilateral has polynomial order of either p = 4 
or p = 6. Table 2 compares the current DG results for the total drag coefficient, C d , and 
separation bubble length, s/d, with prior high-order finite-difference computations and 
experiments. The DG results for both p = 4 and p = 6 are nearly identical, indicating that 
these quantities are converged. The DG results are within about 7% of the experimental 
results, which is a negligible difference given the difficulty of performing measurments 
at such low Reynolds numbers. Comparing the DG results to the prior finite-difference 
calculations of Visbal (1986) yields a difference of about 6% in C d and less than 3% in s/ d. 
Interestingly, Morgan et al. (2002) recently performed simulations using a block-parallel 
version of the same solver used by Visbal (1986) and they report up to 3% difference in 
both s/d sind C d . While the source of the discrepencies between these three codes is not 
known, the DG results are converged with regard to both domain size and resolution. 

5.2. Vortex shedding 

Next, consider unsteady vortex shedding from a circular cylinder. The Reynolds number 
based on diameter and freestream velocity is Re = 100, the freestream Mach number is 
Moo = 0.2 and an isothermal condition is enforced at the cylinder surface with T w = Too- 
We have performed simulations over a range of domain sizes and have investigated both 
h and p-refinement to establish the convergence properties of the method. For brevity, 
we show results only for a relatively large rectangular domain, of size xi € [—15, 30] and 
X2 6 [_30, 30], using 812 quadrilateral elements with a tensor-product basis of Legendre 
polynomials on each element, where the polynomial order varies from p = 5 to 8. We note 
in passing that this domain was found to be sufficiently large to prevent far-field boundary 
influence on the solution. Table 3 documents the convergence of the Strouhal number 
St, peak-to-peak lift coefficient AC;, and average drag coefficient C d with polynomial 
order. We see that even with p- 4 all quantities sire accurate to three significant figures. 
When p = 8 the average drag coefficient is converged to at least 5 significant figures. 
The converged Strouhal number is St — 0.1653 which is in excellent agreement with 
the experimental value of 0.165 (Williamson 1989). For both the steady and unsteady 
cylinder flows, the weak implementation of wall boundary conditions is found to provide 
excellent results, even for rather coarse discretizations. 
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P St 

4 0.1652 

5 0.1652 

6 0.1653 

7 0.1653 

8 0.1653 
Expf 0.165 

Table 3. Convergence with polynomial order for 
fExperimental data is 


AC* C d 
0.6951 1.4104 
0.6953 1.4105 
0.6958 1.4106 
0.6960 1.4107 
0.6960 1.4107 


k shedding from a circular cylinder at Re = 100. 
Williamson (1989) 



Figure 4. Cross-stream (z-y) quadrature grid for the stretched mesh with p = 5, 4, 3. 


6. Fully-developed channel flow 

We now consider fully-developed turbulent flow in a plane channel with coordinates 
x = xi in the streamwise direction, y — x 2 in the wall-normal direction, and z = x 3 in the 
spanwise direction. The flow is assumed to be periodic in the streamwise and spanwise 
directions where the box size is selected so that the turbulence is adequately decorrelated 
in both directions. Coleman et al (1995) provide excellent documentation of DNS results 
for compressible channel flows at low Re r . 

As a first step towards utilizing DG for turbulent flows, we have performed DNS at 
Re r = 100 with a centerline Mach number of M c = 0.3 so that comparisons can be made 
directly to prior incompressible results (see e.g. Kim et al (1987); Moser et al (1999)). 
Following Coleman et al (1995), we use a cold, isothermal wall so that internal energy 
created by molecular dissipation is removed from the domain via heat transfer across 
the walls, allowing a statistically steady state to be achieved. The bulk mass flow is held 
constant by the addition of an xi-momentum source on the right-hand side of (2.1a). 

The computational domain is (4 tt, 2, 47r/3) and this is discretized with an array of 
8x8x8 elements yielding a total of 512 elements. Exploiting the flexibility of the DG 
discretization, we use both h and p refinement to more efficiently resolve flow features 
near the wall. In particular, two wall-normal distributions of elements are investigated. 
We first use a stretched mesh in the wall-normal direction where the grid points are given 
by 


_ tanh(c g (2 j/iVy - 1)) 
i tanh c 8 


3 = 0,1,..., JVy 


(6.1) 


where N y = 8. For this mesh, we reduce the polynomial order away from the wall, 
starting with two layers of p = 5 elements, then a layer of p = 4 followed by a layer 
of p = 3 elements near the centerline. Thus, moving from the bottom wall to the top 
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Figure 5. Discontinuities in instantaneous and averaged mean-flow profiles, Re T — 100: 
(a) instantaneous u, (b) instantaneous p, (c) average u, (d) average p. 


wall, the element order varies like: {5, 5, 4, 3, 3, 4, 5, 5} resulting in a total of 79,488 
degrees of freedom. Note that the flexibility of the DG method makes it possible to 
coarsen simultaneously in all three coordinate directions as one moves away from the 
wall. The crossflow quadrature grid for the stretched mesh is shown in figure 4. We also 
have performed simulations using a uniform h mesh in the wall-normal direction but 
again with variable p order. For this mesh, two p distributions were considered: a low- 
resolution case with p ={5, 5, 4, 3, 3, 4, 5, 5} yielding 79,488 degrees of freedom and a 
high-resolution case with p ={6, 6, 5, 4, 4, 5, 6, 6} resulting in 131,456 degrees of freedom. 
In all cases, we use third-order TVD-RK time advancement with At = 0.0001. This time 
step is a factor of 10 smaller than that typically used in our incompressible code (Collis 
et al. 2000) because the incompressible code treats wall-normal viscous terms implicitly. 
We are currently enhancing our DG code to support implicit time-advancement. 

We also note that computing turbulence statistics from a DG solution requires a sub- 
stantial coding effort, so that currently we compute only mean and rms quantities. Higher- 
order statistics and spectra will be presented in the future. 

We begin by plotting typical instantaneous and average u and p profiles for the 
stretched mesh solution in figure 5. In plotting all the results shown in this paper, no 
smoothing or other postprocessing has been done to remove the discontinuities inherent 
in a DG discretization. Thus, we can clearly see discontinuities in the instantaneous pro- 
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Figure 6. Mean and mns velocity profiles in wall units for the stretched mesh: DG, 

incompressible DNS, law of the wall. 


files, especially in p near the center of the channel. However, after averaging, both the 
stream wise velocity and density profiles are smooth. One of the nice features of DG is 
that if the solution is known to be smooth, then visible jumps in the solution are indica- 
tive of low resolution. Thus, with the stretched mesh, the instantaneous turbulent flow 
near the center of the channel is only marginally resolved, although near the walls even 
the instantaneous profiles appear smooth, indicating good resolution there. However, it is 
important to note that even though the resolution near the centerline may be marginal, 
the mean flow is well represented. 

Evidence to support this claim is given in figure 6 which shows the mean and rms 
velocity profiles in wall units, compared to a reference incompressible DNS at the same 
Re T (Chang 2000). Both the mean and rms velocities are in excellent agreement with 
the incompressible DNS. Likewise, no discernible discontinuities axe observed in either 
the mean or the rms profiles. We recall that the DG discretization uses 79,488 degrees of 
freedom and is formally between 4th- and 6th-order accurate, depending on the local poly- 
nomial order. For comparison, the incompressible DNS uses a hybrid Fourier-Galerkin 
method in the planes and second-order finite-volume method in the wall-normal direction 
and uses 336,960 degrees of freedom after dealiasing. Thus, the DG solution uses a factor 
of 4.2 less degrees of freedom (1.9 if dealiasing is not used in the incompressible case). 

On the stretched mesh, the average slip in the streamwise velocity at the wall is 0.002% 
of the centerline velocity where the first collocation point is A y+ = 0.7 from the wall, f 
To determine how the weak wall boundary condition influences the solution at coarser 
resolution (near the wall) we now consider results using a uniform mesh in the wall- 
normal direction. Figure 7 shows the mean streamwise velocity profiles in wall units as 
compared to the reference incompressible DNS, for both the low- and high-resolution 
cases. Interestingly, we see that the profiles are in excellent agreement with the reference 
solution. Such overlap clearly indicates that the mean shear stress is well predicted in 
both cases. However, careful examination of figure 7 does show that the law of the wall 
u+ = y+ is not perfectly satisfied at small y + because the flow slips at the wall. For the 
low resolution case, the slip velocity is 1% of the centerline velocity with A y+ = 2 while 
for the higher resolution case there is 0.68% slip with A = 1.6. As expected, as near- 
wall resolution is increased, the amount of slip is reduced as the enforcement of the wall 

f For reference, the centerline velocity is approximately 16u r at Re r = 100. 
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Figure 7. Effect of slip on mean velocity profile for the uniform mesh at two resolutions: DG, 

incompressible DNS, law of the wall. 



Figure 8. Effect of slip on rms velocity profiles for the uniform mesh at two resolutions: 

incompressible DNS. 


DG, 


boundary condition improves (this is especially true for the stretched mesh) . Importantly, 
the mean shear and the majority of the mean velocity profile are well predicted even for 
the lowest-resolution case when A = 2, which is less than many RANS models allow. 

Similar behavior is found for the rms velocities, as shown in figure 8 for the low- and 
high-resolution uniform-mesh cases. One can clearly see the slip in the streamwise rms 
velocities at the wall. For the low-resolution case u+ m$ = 0.65 at the wall, while for the 
high-resolution case u+ ms = 0.48 at the wall. For reference, the stretched-mesh solution 
has u+ ms = 0.0062 at wall. Again, as resolution is increased in the near-wall region, the 
no-slip boundary condition is enforced to a higher accuracy. Importantly, with the ex- 
ception of a region very close to the wall, both the mean and rms profiles throughout 
the channel are well predicted for all cases. Our prior experience with hard boundary 
conditions has shown that mean shear and rms profiles (i.e. turbulence production) are 
incorrectly predicted at low resolutions. Conversely, by enforcing the wall boundary con- 
ditions weakly through the numerical fluxes, the influence of the wall on the flow is 
correctly simulated in the form of net shear stress and turbulence production, even at 
resolutions for which the wall boundary values are inaccurate. 
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7. Conclusions 

A discontinuous Galerkin method is formulated and implemented for simulation of com- 
plex, turbulent, compressible flows. The implementation is validated for both steady and 
unsteady separated flow over a circular cylinder, with results in excellent agreement with 
prior computations and/or experiments. An important feature of discontinuous Galerkin 
is the ability to enforce Dirichlet boundary conditions weakly, through numerical fluxes 
at the wall. The advantages of this approach are demonstrated for a simple advection- 
diffusion problem, where it it shown that enforcement of a weak boundary condition leads 
to a significant reduction in oscillations in the computed solution, resulting in a factor 
of 4 times less error in the L 2 norm. Applying DG to simulate fully-developed turbulent 
flow in a plane channel at low Reynolds number Re r = 100 leads to results in excel- 
lent agreement with a reference incompressible DNS. The advantage of weak Dirichlet 
boundary enforcement is also demonstrated for this flow, where it is shown that accu- 
rate profiles of net shear stress, as well as mean and rms velocity, are obtained at low 
resolution — even resolution for which there is significant slip at the wall. In this context, 
weakly enforced wall boundary conditions may play a useful role in wall modeling for 
large-eddy simulation, where the wall-model is given by a particular numerical flux used 
at the wadi. 
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